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ABSTRACT 

We investigate the rotational evolution of solar-like stars with a focus on the inter- 
nal angular momentum transport processes. The double zone model, in which the 
star's radiative core and convective envelope are assumed to rotate as solid bodies, 
is used to test simple relationships between the core-envelope coupling timescale, r c , 
and rotational properties, like the envelope angular velocity or the differential rotation 
at the core-envelope interface. The trial relationships are tested by fitting the model 
parameters to available observations via a Monte Carlo Markov Chain method. The 
synthetic distributions are tested for compatibility with their observational counter- 
parts by means of the standard Kolmogorov-Smirnov (KS) test. 

A power-law dependence of r c on the inner differential rotation leads to a more 
satisfactory agreement with observations than a two- value prescription for r c , which 
would imply a dichotomy between the initially slow (P rot > 3 d) and fast (P ro t < 3 d) 
rotators. However, we find it impossible to reconcile the high fraction of fast rotators 
in a Per with the rotation period distributions in stellar systems at earlier and later 
evolutionary stages. This could be explained by local environmental effects (e.g. early 
removal of circumstellar discs due to UV radiation and winds from nearby high-mass 
stars) or by observational biases. 

The low KS probability that the synthetic and observed distributions are not 
incompatible, found in some cases, may be due to over-simplified assumptions of the 
double zone model, but the large relative uncertainties in the age determination of very 
young clusters and associations are expected to play a relevant role. Other possible 
limitations and uncertainties are discussed. 

Key words: Methods: numerical - stars: rotation - stars: interiors - stars: late-type 
- stars: pre-main sequence 



1 INTRODUCTION 

Rotation is an important parameter for solar-like stars. 
It is an essential ingredient of th e dynamo process (e.g. 
iBrandenburg fc Subramanianll2~005t ). which in turn powers 
magnetic activity, it can influence the surface abunda nce of 
light elements through rotationally-induced mixing ( Zahn 
1 19931 ) . and it is related to the formation and evolution o f 
planetary systems (see, e.g.. lChang et al ] |2010l : lLanzall20ld ). 

Rotation period measurements in open clusters and 
young associations provide a valuable source of information 
on angular momentum evolution of late-type stars. In recent 
years, photometric measurements of the stellar rotation pe- 
riod, based on the rotational modulation of the light curve 
induced by photospheric starspots, have replaced the mea- 
surements based on the v sini as derived from spectral line 



profiles. This has considerably improved both the number 
and the precision of available rotation periods. Such obser- 
vational constraint, however, has not been exploited in its 
full potential yet. In fact, most of the analyses carried out so 
far have considered only limiting cases such as the subsets 
of the faster and slower rotators at given ages. 

During the pre-main sequence (PMS) phase, solar- 
mass stars undergo a global contraction. However, for the 
first few Myr, the interaction with a circumstellar disc 
drains angular momentum from the star, thus delaying 
its spin up for the (variable) duration of the disc life- 
time . _Thi5_rjrocess_J^_j2^^ 



Collier Cameron fc Campbell 1 19931 : IShu et alj Il994l : 



iMatt et all l2010h and is usually modelled by means of 
very simplifie d assumptions, e.g., the disc- locking hypothesis 
l|Konig]||l99ll) . Once on the main sequence (MS), the rota- 



2 F. Spada et al. 



tional evolution is driven by the loss of angular momentum 
through a magnetised stellar wind (Weber fe DavTsl 1 19671 ; 
iKawalerl 119881 ; lchabover et~aj]ll995al lbh. 

Surface period measurements are also an indirect 
probe of the internal rotation profile of these stars. Phe- 
nomenological modelling of MS angular momentum evo- 
lution ( the so-called double zone m o del, DZM in the fol- 
lowing: iMacGregor fc Brennerl Il99ll ; iKeppens et~al] 1 19951 ; 
lAllairJll998l : lBouvieril2008l . see also Sect.GZE}, assuming that 
the radiative core and the convective envelope rotate rigidly 
at all ages, has successfully reproduced the observations 
provided a certain amount of differential rotation between 
the core and the envelope is allowed. This is parameterised 
through a coupling timescale r c , which determines the rate 
of angular momentum exchange between the two regions. 
From the comparison of synthetic rotational tracks and the 
upper and low e r perc entiles of the observed period distribu- 
tions, iBouvieri {2008) concluded that r c is of the order of 10 
Myr for the stars that begun their evolution on the ZAMS as 
fast rotators (i.e., with an initial rotation period of ~ Id), 
while it is remarkably longer (t c ~ 100 Myr) for initially 
slow rotators (i.e., with initial periods ~ 10 d). A possi- 
ble connection between this dichotomy and the presence or 
absence of a planetary system orbit ing the star is further 
discussed in lEggenberger et all |2010| . An understanding of 
the physical processes that give rise to such a dichotomy 
and eventually establish a nearly solid-body rotation within 
the age of the Sun (as obse rved through helioseismology, see 
e.g. iThompson et alj|2003h is still lacking. 

Our aim is twofold: on the one hand, we try to extract 
the whole statistical information from the observed rotation 
period distributions. In recent years, rotation period mea- 
surements have improved significantly both in number and 
in precision. In spite of the age spread within cluster mem- 
bers, uncertainties on cluster age estimates and other sys- 
tematic effects, good sample completeness and homogeneity 
has already been achieved. We exploit this information com- 
paring the observed period distributions with synthetic ones, 
constructed by evolving the Orion Nebula Cluster (hereafter 
ONC) distribution by means of the DZM. The modelling of 
the period distributions based on the DZM requires some 
prescription on the dependence of r e on the an gular veloc- 
ity, to comply with the findings of lBouvi er (2008) mentioned 
above. On the basis of intuitive physical arguments, we for- 
mulate some hypotheses on the scaling of r c . The set of pa- 
rameters in the DZM that produces the best fit of the syn- 
thetic distributions to the observed ones is determined by 
means of an iterative procedure based on a Markov Chain 
Monte Carlo (MCMC) sampling. This suffices to the pur- 
pose of comparing the different hypotheses on r c and to test 
the overall capability of the DZM to satisfy the latest obser- 
vational constraints. 

A direct comparison of theoretical modelling of stel- 
lar rotation with observed period dist ributions has alread y 
been attempted in the past by, e.g., iTinker et af] (2002). 
They evolved the ONC period distribution to the age of the 
Pleiades and the Hyades by means of the Yale Rotating Evo- 
lution Code |Guenther et al.lll992j ). mainly aiming at deter- 
mining the parameters in the wind braking law. The present 
study is novel in many respects, because we rely upon a sim- 
pler modelling, focussing on the coupling timescale. More- 
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Figure 1. Observed period distribution for the 57 stars of ONC 
in the selected mass range. 



over, a much larger sample of observed period distributions 
is used. 



2 ROTATIONAL PERIOD DATA 

We use the rotation period distributions extracted from pho- 
tometric surveys of stellar open clusters and young associ- 
ations, as available in the literature (see Table [1] for the 
references). 

The rotational data for the associations TW Hydrae, 
/3 Pictoris, Tucana-Horologium, Columba, Carina (age < 
30 Myr), AB Doradus (age ~ 120 Myr) are taken from 
iMessina et"afl (l2010al) and for th e cluster M 11 (age ~ 230 
Myr) from lMessina et al.l l|2010bT ) and are used for the first 
time for this kind of investigation. Remarkably, these young 
associations fill an age gap previously present between Tau- 
Aur (~ 6 Myr old) and a Per (~ 70 Myr) in the mass range 
0.9 < M,/M Q < 1.1. 

The ages reported here are the most recent literature es- 
timates. In particular, the association AB Dor is considered 
al most coeval with the Pleiades, according to the analysis 
of lMessina et al l (|2010bh . who revised a previously proposed 
age of 70 Myr. A remarkably extensi ve rotational dataset fo r 
the Pleiades have been produced bv lHartman et all (|2010l ). 

The rotational evolution significantly depends on stel- 
lar mass, owing to different durations of the PMS phase, 
convective envelope depths, etc. As we are primarily in- 
terested in the rotational evolution of solar analogues, a 
subsample of each period distribution in the mass range 
0.9 ^ M*/Mq ^ 1.1 was selected. For very young clus- 
ters (i.e. age ^ 10 Myr), which are affected by a substan- 
tial reddening, the effective temperature of the stars was 
determined on the basis of their spectral type as given by 
iRebull et all (120031 ). For MS objects, the dereddened (B-V) 
and (V — /) colours were employed for the same purpose. In 
both cases, the determination of the mass from T e g or the 
colour index was perfor med by direct compa rison with the 
theoretical isochrones of iBaraffe et al.l (|l998l ). 

The ONC period distribution, providing the initial con- 
ditions for our rotational modelling, is shown in Figure [1] 
In the restricted mass range which is of interest to us, the 
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Table 1. Names, estimated ages and number of stars in the 
mass range 0.9 ^ M*/Mq ^ 1.1 for the stellar systems used in 
this work. Stellar associations arc in bold. 



Object name 


Age [Myr] 


# stars 


Reference 


ONC 


2 


57 


[a] 


NGC 2264 


4 


52 


[a] 


NGC 2362 


5 


43 




Tau-Aur 


6 


25 




TW Hya, /3 Pic 


10 


12 




Tuc-Hor, Col, Car 


30 


46 


© 


a Per 


70 


43 


El 


AB Dor, Pleiades 


120 


109 


|c]|e] 


M 50 


130 


62 


13 


M 35 


150 


88 




M 34 


220 


17 


1 


M 11 


230 


15 


m 


M 37 


550 


89 





(a) : compilation, see iRebull et al . (2003) and references 
therein; 

llrwin et all feOOSl) ; 
iMessina etafl [|2010ah ; 
co mpilation, see Messina et al. | fcOOlf) , | Messina et all 
and references therein; 



(b) 
(c) 
(d) 



(e) 
(f) 
(g) 
(h) 

(i) 

Ci) 



lHartman et al.l j2010h : 
llrwin et all j200Sh ; 
iMeibom etaD (2009); 
llrwin et all l|2006l) : 
IMessina et~all ll2010bl) ; 
lHartman et al.l ||2009|) . 



distribution is unimodal, with a fairly pronounced peak of 
fast rotators (i.e., with P rot ^ 4 days). 



3 NUMERICAL PROCEDURE 

We construct a series of synthetic rotation period distribu- 
tions at the ages reported in Table [T] Each period in the 
ONC solar analogue subsample is evolved forwards in time 
with the DZM, calculated for a 1 Mq star. The synthetic 
distributions are tested for compatibility with their observa- 
tional counterparts by means of the sta ndard Kolmogorov - 
Smirnov (KS) test (see, e.g., sect. 14.3 of lPress et al.l ll992). 
The input to our numerical procedure is an initial guess 
for the parameters appearing in the equations |l} of the 
DZM. The parameter space is then sampled by means of the 
MCMC method, which performs random jumps in each pa- 
rameter, guided by a probabilistic acceptance rule (see Sect. 
13. 3|) . The best-fitting set of model parameters is obtained in 
output. 



3.1 Parameters in the DZM 

Our simple model for the rotational evolution of a solar-like 
star assumes that the radiative core and the convective en- 
velope both rotate as sol i d bodies (|MacGregor &; Brenner! 
Il99ll ; iKeppens et"ai1ll995l ; lAUairJll998h . The envelope is ex- 
pected to be strongly coupled by the very efficient angular 
momentum redistribution associated with turbulent convec- 
tion. On the other hand, a large scale magnetic field, even 



as weak as < 1 G, can maint ain a condition of r igid rotation 
in most of the core (see, e.g. lSpada et~afll2010h . 

At any given time t, the angular momentum of the core, 
J c , and of the envelope, J e , are: 

J c (t) = I c (t) Q c (t) ; J e (t) = hit) « e (t), 

where I c> I e and f2 c , Q e are the moments of inertia and the 
angular velocities of the core and the envelope, respectively. 
Our evolution begins at the ONC age, i.e. to ~ 2 Myr. Stars 
are fully convective until to, thus justifying a solid-body ro- 
tation as our initial condition, i.e. fl e (to) = Q, c (to)- 

The evolution of the angular momenta is governed by 
two coupled differential equations: 



dJ c 
~~dt 

dj e 
~~dt 



AJ 



2 2 dM c 

3 c dt 



= +- 



AJ 



2 2 

3 Rc 



dM c 
dt 



O,. 



(1) 



dJ e 

~dt 



where R c and M c are the radius and mass of the radia- 
tive core. The stellar structure quantities and their time 
deriva tives are taken from PMS and MS evolutionary mod- 
els by iBaraffe et al.l (| 19981 ') . In equations ([1]) , the core and 
the envelope exchange an amount of angular momentum 

AJ = e r c — - at a rate determined by the coupling 

Ie+Lc 

timescale r c . During the PMS, the growth of the core at the 
expenses of the envelope produces an angular momentum 
transfer, accounted for by the second term in the right-hand 
side of both equations (JTJ. Because of the magnetised wind, 

the envelope loses angular momentum at a rate — ^ 

dt 



for w hich the following parametric formula is used (|Kawalerl 
GUI): 



dj e 

~~dt 



= K U 



R*/R@ 
M*/Mz 



1/2 



min(ntnL (2) 



where K w and f2 sa t determine the braking intensity and sat- 
uration threshold, respectively, while R*/Rq and M»/M© 
are the stellar radius and mass in solar units. Finally, the 
angular momentum loss due to the interaction with the disc 
during early PMS evolu tion is taken into account with the 
disc-locking hypothesis (|Konigllll99ll ). i.e.: 



dQ e 
~~dt 



while t ^ Tdi s 



In other words, the net effect of the interaction with the 
disc is that of keeping the surface angular velocity of the 
star constant for the whole disc lifetime. 

The solutions of equations |TJ) depend on the parame- 
ters K w , S7 sat , t c and Tdi sc - In applying the DZM to a whole 
distribution of rotators, the wind law parameters K m and 
f2 sat can be assumed to be the same for all the stars; for 
t c and Tdisc, on the contrary, such an assumption is rather 
crude. It has already been established that it is impossible 
to reconcile the lower and upper envelopes of the rotation 
perio d distributions with the same value of t c (e.g. iBouvierl 
2008). This is dealt with by simple assumptions on the de- 
pendence of t c on stellar rotation, as discussed below (see 
Sect. 13.2] ). A realistic treatment of the star-disc interaction is 
a complex task. Even under the disc-locking assumption, the 
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value of Tdisc for each star in the sample should be extracted 
from a distribution. 

Although information on the distribution of disc life- 
time, as well as on its dependence on the presence of high 
mass stars in t he neighbourhood, is becoming increasingly 
availa ble (e.g. iHernandez et al.l 120081 ; iKennedv fc Kenvonl 
2009), we refrain from exploiting it in full details so to keep 
our model as simple as possible. 



3.2 Hypotheses on the coupling timescale 

We consider the following prescriptions on the dependence 
of t c on stellar rotation: 

• A two-valued function of the form: 



10 Myr 

TQ 



otherwise 



(3) 



This is the si mplest implementation of t he early empiri- 
cal fi nding (e.g. iKrishnamurthi et ailll997l . see also iBouvierl 
2008) that fast rotators are well described by a nearly solid- 
body inner profile, while a sizeable decoupling is necessary 
for slow rotators. The value of 10 Myr used for initially 
fast rotators is so short as t o ensure nearly instanta neous 
coupling at all ages (see also iDenissenkov et al.| [2010). The 
timescale for slow rotators, To, and the initial angular ve- 
locity threshold, fi C rit 5 are treated as free parameters to be 
determined by means of the MCMC procedure. 

• A power law dependence on the instantaneous amount 
of differential rotation: 



r c {t) = T 



AQ(t) 



(4) 



where Afi(t) = Q c (t) — Q e (t) and Af2 is the inner differ- 
ential rotation of the present Sun, assumed as a reference 
value (we used AS1 = 0.2 fi ). 

Note that this choice produces a time-dependent r c . This 
is an attempt to account for the effect of an enhanced an- 
gular momentum transport due to turbulent viscosity. The 
viscosity enhancem ent could b e the result of either a rota- 



tional mixing (e.g. IZahnl Il993h or of hydromagnetic insta- 
bilities (|Spruitlll999l . l2002h ~ powered by the available kinetic 
energy stored in the differential rotation. The form of equa- 
tion Q, with a > 0, is chosen to produce a shorter coupling 
timescale (which mimics a greater effective viscosity, ensur- 
ing a stronger coupling) when a greater differential rotation 
is present, to and a are regarded as free parameters. 

Other simple prescriptions, e.g. a dependence of r c on 
the surface angular velocity £l e {t), failed to give suitable 
results and are not discussed in the following. 

3.3 The MCMC method 

Our application o f the MCMC method (see, e.g. 
iTegmark et al .120041 ) constructs a sequence of values for each 
parameter of the DZM, approaching the best fitting region 
of the parameter space quite rapidly (burn-in phase, usually 
a few hundreds of steps), and then sampling this restricted 
neighbourhood by performing a random walk within it. 

The routine is initialised providing a starting guess for 
the parameters. At the beginning of each step, the current 



Table 2. Best fitting parameters estimated applying the MCMC 
method. The units of K w are g cm 2 s. 





two-valued t c 


power law t c 




(cf. equation [3l 


(cf. equation |4)l 


(K w ) 


= (5.40 ±0.054) ■ 10 4 ' 


(K w ) = (5.97 ±0.13) ■ 10 4 ' 


(to 


) = 128 ± 3.84 Myr 


(t ) = 57.7 ±5.24 Myr 




it) = 3.89 ± 0.14 H 


(a) = 0.076 ± 0.02 



set of parameters is used to evolve forwards in time the 
ONC period distribution to each of the ages in Table [T] by 
means of the DZM. A KS test is then performed comparing 
the synthetic period distributions with their observational 
counterparts. Calling Pks,» the probability of the synthetic 
and observed distributions to be significantly different for 
a certain cluster i and viewing the individual KS tests as 
independent from each other, we define the quantity Q by 
adding these probabilities together: 



KS,i , 



with N being the number of rotation period distributions 
used. Two additional constraints come from the rotational 
properties of the present Sun, namely: 



n c (t Q ) - n e (t ) 
n e (t ) 



<0.2, 



(5) 



which are taken into account by suitable additive terms in 
Q. The MCMC procedure is used to iterate the model pa- 
rameters towards the minimum of Q, which corresponds to 
the best fit. 

After an evaluation of Q, the following step begins with 
a tentative random jump in each parameter. The new set 
of parameters just generated is retained in the chain or 
rejected following the stochastic Met ropolis-Hastings rule 
(|Metropolis et al.il 19551 : lHastingsl Il970h : 

AQ < : accept; 

AQ > : accept with probability oc e~ A ® , 

where AQ is the difference between the current value of Q 
and its value at the last accepted step. This rule ensures a 
quick descent towards the global minimum, with a non-zero 
probability of escaping from possible local minima encoun- 
tered during the path. 



4 RESULTS 

The procedure described above was applied to determine the 
best values of K w , to and f2 C rit, in equation @, or a, in the 
case of equation Q, respectively. The parameters f2 sa t and 
Tdisc were held fixed. We used the following fiducial values: 



1 2^ 



5.5 fir. 



Tdh 



5.8 Myr. 



These choices are motivated by independent observational 
estimates, namely the angular velocity saturation thresh- 
old for the X-ray em ission of late-type stars (see figure 3 of 
iPizzolato et alj|2003h for fi sa t, and the age at which about 
10 % of stars still possess a circumstellar disc (see figure 
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Figure 2. Values of r c determined by the power-law prescription 
with the best-fitting tq and a (see Table [2} for two initial veloc- 
ities representative of slow (dotted line) and fast (dashed line) 
rotators, respectively. 



1 of iMamaiekl [20091 for raise- This is also in Rood quan- 
titative agreement wi th earlier works (e.g. iBouvie rl |200S| ; 
llrwin fc Bouvierll2009l ). 

The best-fitting values of K w , to and Q c nt were calcu- 
lated from the portion of the MCMC chain that effectively 
sampled the region of the parameter space around the min- 
imum (in both cases, we had > 10 4 steps). For each param- 
eter, we used the mean value and the standard deviation 
as estimates of the best fitting value and its uncertainty, 
respectively. The results are summarised in Table [5] Note 
that the value of K w found, of the order of 10 47 gem 2 s, 



is in good agreem ent with previous works (see, e.g. Kawalerl 
ll988l . lAllainlll998l ). 

To allow the comparison of the two different prescrip- 
tion for r c , we show in Figure [2] the evolution of r c in the 
power-law case for two initial values of f2 e representative of 
slow and fast rotators. Remarkably, the best-fitting param- 
eters produce a r c > 40 Myr for fast rotators, significantly 
longer than 10 Myr, which was used in the two-valued case 
to ensure nearly solid-body rotation at all ages, as suggested 
by previous works. 

A comparison of the quality of the fitting in the two 
cases, based on the individual KS tests, is presented in Fig- 
ures [3] and [4] which show the synthetic and observed cu- 
mulative distributions (CDF) for the various clusters. The 
KS compatibility probability, Pks.i, i-e., the probability 
of the two distributions for cluster i of being not incom- 
patible, is reported in each panel along with the number 
of stars in the respective observed distribution (note that 
fks.i = 1 — Pks,i)- The smaller the value of Pks,i, the worse 
the fitting between the calculated and the observed distri- 
butions for cluster i. 

The photometric periods on which our analysis is based 
are Fourier-derived from the rotational light modulation in- 
duced by starspots. The amplitude of such a modulation 
falls below the minimum value currently detectable from the 
ground (~ 0.01 mag) for MS stars with P rot > 10 d (see also 
Sect.JSJ). As a consequence, the photometric rotation period 
of slow rotators is hardly derived and the observed distri- 
butions may suffer a lack of completeness beyond this value 



of the rotation period. To account for this effect, periods 
in the synthetic distribution greater than the maximum of 
the corresponding observed distribution are not considered 
when calculating Pksa and Pks.i- 

A good agreement (Pks,i > 0.37) is achieved for the 
first ~ 10 Myr in both Figures. The first two panels refer to 
clusters younger than the assumed Tdi sc ; as a consequence, 
the corresponding synthetic distributions show no evolution 
with respect to ONC, by virtue of the disc-locking con- 
straint. For ages greater than raise but lower than ~ 30 Myr, 
the spin-up dominates due to the stellar contraction. At 
6 and 10 Myr, therefore, the agreement with observations 
mainly rely on Tdi sc and on the stellar structure and evo- 
lution model, while it is quite insensitive to r c and K w . 
For ages older than approximately 30 Myr, the agreement is 
strongly dependent on K w and on the model assumed for t c . 
For the two- value r c model, the only satisfactory agreement 
is found for M 34 (220 Myr). For the power-law r c model, a 
satisfactory agreement is found for the Pleiades (120 Myr) 
and M 37 (550 Myr). We note, however, that the CDFs de- 
rived from the observations do not follow an evolution as 
regular as expected from the models. The a Persei CDF, in 
particular, displays a sharp rise at short rotation periods due 
to an excess of fast and ultra-fast rotators. Such a sharp rise 
is not seen in any other CDFs, except, perhaps, in M 34 and 
M 11, although the number of periods available is too low 
to draw any definitive conclusions. The two-value r c model 
can account for this sharp rise, but then fails to reproduce 
the smooth CDF shape of AB Dor/Pleiades, M 50, etc. The 
power-law r c model, on the other hand, cannot reproduce 
any sharp rise in the CDF and therefore fails to reproduce 
the a Persei, M 34 and M 11 CDFs. 

Without any angular momentum transfer from the core 
to the envelope, the maximum angular velocity is expected 
at about the Tuc-Hor / Col / Car age, i.e., « 30 Myr (see 
also Figures 0and[5]). However, this is not observed because 
the evolution from the Tuc-Hor / Col / Car age (30 Myr) to 
a Persei (70 Myr) appears consistent with a continuation of 
the spin-up after the ZAMS. The spin-up of fast and ultra- 
fast rotators from the Tuc-Hor / Col / Car age (30 Myr) to 
the a Persei (70 Myr) age seems even faster than for slower 
rotators. In the DZM framework, such a behaviour cannot 
be reproduced unless, for the fast and ultra-fast rotators, 
a remarkable amount of angular momentum is transferred 
from the core to the envelope within a very short timescale 
(~ 30 — 40 Myr). In this case, however, very little angular 
momentum would be left in the core for the later evolution. 

After the peak in the angular velocity, a generally mono- 
tonic spin-down is expected. The fast rotator excess in M 34 
(220 Myr) and Mil (230 Myr), however, hampers the pos- 
sibility of obtaining a satisfactory fit at all ages. In fact, the 
percentage of observed fast rotators (P ro t ~ 1 - 2 d) tends 
to decrease from a Per to M 35, but then increases abruptly 
in M 34 and M 11. Such a distribution evolution is clearly 
inconsistent with the DZM alone. 

Despite such limitations, which will be discussed in 
some more details in Sect.0 the power-law r c model pro- 
duces a more satisfactory fit to the data, particularly since 
it reproduces the shape of most of the observed CDFs also 
giving a high Pks.i value for AB Dor / Pleaides and M 37. It 
requires, however, an alternative explanation for the excess 
of fast and ultra- fast rotators in a Per, M 34, and Mil. 
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Figure 3. Synthetic (grey line) and observed (black line) rotation period cumulative distributions for the two-valued r c case. The 
complementary probability of being drawn from significantly different distributions, as calculated by means of the KS test, is reported 
in each panel together with the number of objects in the observed distribution. 
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Figure 4. Synthetic (grey line) and observed (black line) rotation period cumulative distributions for the power law r c case. The 
complementary probability of being drawn from significantly different distribution, as calculated by means of the KS test, is reported in 
each panel together with the number of objects in the observed distribution. 
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Two rotational tracks, with initial conditions represen- 
tative of slow (fi e (io) = fic(to) = 3 Qq) and fast rotators 
(fi e (to) = O c (to) = 25 n ), as calculated with the DZM 
and prescriptions |3} or J|[ for r c , are compared with the 
observed distributions and the solar-age angular velocity in 
Figures and [5J respectively. Assigning a constant value 
to t c for the whole rotational evolution results in slow and 
fast tracks that intersect with each other at ages later than 
about 1 Gyr (see Figure [5]). This is another consequence of 
the high percentage of fast rotators in a Per, which imposes 
a fast angular momentum transfer from the core to the enve- 
lope just after the ZAMS. This diminishes dramatically the 
core angular momentum reservoir available for later stages, 
producing a faster spin-down for fast rotators. On the other 
hand, the power-law prescription for r c leads to a remark- 
ably greater internal differential rotation for the fast rotators 
and a somewhat lower internal differential rotation for the 
slow rotators than the two- value prescription. 

The solar constraint (equation [5} is reasonably matched 
by the power-law prescription. The large differential rotation 
developed by fast rotators is due to the relatively large value 
of t c , as previously noted. In contrast, the poor matching 
found for the slow rotators in the two-value case may seem 
to be in contradiction with previous results (e.g. iBouvierl 
2008). Note, however, that in the present work an attempt is 
made to fit the rotation period distributions, not just their 
percentiles. A value of r c ~ 130 Myr, which is the result 
of our best-fitting procedure (and as such is the outcome 
of taking into account all the period distributions plus the 
solar constraint), turns out to be too large to recover the 
almost rigid rotation at the age of the Sun (cf. also Figure 
[2] for the power-law case). The slightly lower value of K w 
with respect to the power-law case (see Table [2]), moreover, 
explains the excess of rotation at t@. Indeed, this disagree- 
ment is partly relieved if the a Per distribution, which is 
likely prone to observational biases (see Sect.[5|, is excluded 
from the fitting procedure. A run performed for compari- 
son purposes without the a Per constraint in the two-value 
case gives KS probabilities above 20 % for the intermediate 
age clusters (i.e. M 50 to M 11) and different values of the 
best fitting parameters, namely K w ~ 6.3 • 10 47 g cm 2 s and 
t c ~ 80 Myr. A similar test for the power-law case results in 
negligible differences both in KS probabilities and in best-fit 
parameters, consistently with a moderate sensitivity to the 
a Per constraint. 

In conclusion, the power-law prescription is more suited 
to fit the shape of the rotational distributions, while the two- 
value models better reproduces the basic features of the two 
extreme populations (i.e. the slow and fast rotators). The 
two prescriptions also show different sensitivity to peculiar 
characteristics (either real or due to biases) of the distribu- 
tions. Clearly, none of the two models is able to completely 
describe rotational evolution. 



5 DISCUSSION 

We modelled the evolution of rotation period distributions 
from the PMS to the solar age, using the available observa- 
tions to determine the best-fitting parameters of a modified 
version of the DZM. We aimed at including in the analysis 
the whole statistical information contained in the observed 



distributions as available so far; at the same time, we tested 
the possibility of including in the DZM some kind of depen- 
dence of the coupling timescale on the rotational state of 
the star. A numerical procedure was implemented to evolve 
forwards in time the solar analogues' subsample of the ONC 
period distribution, used as our initial condition. The agree- 
ment of the synthetic distributions with the observed ones 
was evaluated by means of the KS test. This procedure was 
used to select the values of model parameters ensuring the 
best fitting of the observational constraints. 

A prescription for the dependence of r c on the rota- 
tion of the star is required to fit the evolution of the whole 
period distributions. It is already known, in fact, that the 
slow and fast percentiles of the period d istributions req uire 
significantly different values of r c (e.g. IBouvierl 120081 ). In 
formulating these prescriptions (see equations [3] and [4} we 
relied on basic physical intuition, keeping the number of free 
parameters to a minimum. Equation ((3} is the simplest way 
to mimic the dichotomy in r c between slow and fast rotators 
as suggested by earlier works. It implies that different t c val- 
ues are set at early stages by some mechanism that quickly 
distinguishes between two separate rotation regimes, but is 
poorly sensitive to internal differential rotation. A different 
magnetic field configuration withi n the stellar interior might 
accou nt for such a behaviour (see ISpru it 1999; Spad a~t al.l 
120101) . Equation (0, on the other hand, assumes that the 
available content of rotational kinetic energy, as measured in 
the DZM by the amount of differential rotation, can trigger 
angular momentum transport processes with an efficiency 
parameterised by the exponent a. Other formulations, with 
a power law dependence of r c on the surface angular velocity 
Q, e (t), with or without saturation mechanisms like that in 
equation ((2|), proved to be completely unsuccessful. 

A major limitation of this work comes from the uncer- 
tainty in the ages of the clusters (e.g. lMavne fc Navlor|[2008l . 
and references therein). This is particularly relevant for the 
very young clusters and for the stellar associations, which 
are also more likely to be prone to remarkable internal age 
scattering. 

During the PMS, the fitting is dependent essentially on 
raise and on the radius contraction as predicted by the evolu- 
tionary model. Considering the uncertainties in age and the 
crude approximation on the role of the circumstellar disc, 
the fitting can be considered quite satisfactory in this part 
of the evolution. After the ZAMS, when the rotational evo- 
lution becomes very sensitive to angular momentum transfer 
from the core to the envelope and to wind-braking, the ob- 
served CDFs do not show a regular evolution as expected 
from the DZM assumptions. The a Per CDF, in particular, 
shows a marked (ultra) fast-rotators excess, which is diffi- 
cult to reconcile with the CDF of younger systems, either 
assuming a two-value or a power-law prescription for r c in 
the DZM. This distribution feature is not seen in any other 
cluster or association, except perhaps M 34 and M 11. The 
consequences on the model fitting have been discussed in 
Sect.EI 

Obviously, the reasons for such a discrepancy could be: 
a) the DZM is inadequate or contains oversimplified as- 
sumptions; b) the stellar rotation period distributions in 
open clusters and young associations may be affected by 
characteristics of the (parent) cluster as a whole, particu- 
larly those leading to the evaporation of circumstellar discs 
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Figure 5. Evolution of two rotators, according to the DZM as modified by the two-valued t c prescription in equation ((3j and using the 
parameters on the left of Table [2] with initial conditions representative of the fast and slow groups (f! e (t) are plotted as solid lines and 
Q c (t) as dash-dotted lines). Observed rotation period distributions (diamonds) and the Sun (triangle) are also plotted for comparison. 




Figure 6. Evolution of two rotators, according to the DZM as modified by the power law r c prescription in equation Q and using the 
parameters on the right of Tabled with initial conditions representative of the fast and slow groups (Q e (t) are plotted as solid lines and 
£1 c (t) as dash-dotted lines). Observed rotation period distributions (diamonds) and the Sun (triangle) are also plotted for comparison. 
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by U V radiation and wind from neighbour high-m ass stars 
(e.g. iHernandez et all l200St iGuarcello et al.ll2010l ); c) the 
observed period distributions are incomplete and/or biased. 

Indeed, our model does not include differences in the 
rotational history of clusters' stars that may arise from the 
presence of high mass stars, multiple or prolonged star for- 
mation events, the compactness and richness of the clusters, 
and the stars' dispersal into the field. Therefore, despite the 
DZM may contain the essential physics to describe the stel- 
lar angular momentum evolution, a satisfactory modelling 
may require taking into account the peculiarities of the clus- 
ter to which the stars belong. In the case of a Per, for in- 
stance, high mass stars might have influenced the rotational 
history of a fraction of low mass stars by evaporating their 
circumstellar discs at an early stage. Without an efficient 
disc-locking, these stars may have reached the ZAMS (first 
30 Myr) with a very high rotation and have had no time to 
spin-down significantly in the subsequent ~ 40 Myr. 

Other possible sources of uncertainties in the period 
distribution include clusters contamination by field stars, 
the presence of close binaries whose high rotation rate is 
maintained by the tidal synchronisation between the com- 
ponents, and observational biases. Concerning these latter, 
photometric monitoring with a limited temporal extension 
(generally of 10-15 d as in the cases of M 11, M 34, M 35, 
M 50, NGC 2362 and of numerous members of a Per) pre- 
vents us from detecting the slow rotation periods, making 
the derived distributions incomplete. Moreover, for MS stars 
with P ro t > 10 d the amplitude of the optical light curve is 
very small (< 0.001 — 0.01 mag), so their rotational modula- 
tion cannot be measured by ground-based ph otometry (see, 
e.g., iRadick et a"Dll99Sl ; iMessina et al.ll2003l ). The absence 
of rotators with periods longer than ~ 10 d in the observed 
distributions shown in Figures [3] and [4] can therefore be due 
to such limitations, although this is unlikely to hamper the 
analysis to a significant extent as the percentage of rotators 
with P ro t > 10 is predicted to be rather low at young ages. 
On the other hand, a poor sampling from ground-based ob- 
servations (for example, 1-2 d in the case of young associa- 
tions considered in our analysis) makes it hard to detect the 
ultra-fast rotators, making the derived period distributions 
incomplete also on the fast side. Despite high-precision CCD 
photometry and multi-site monitoring campaigns may over- 
come most of these difficulties, the data available to date 
are still rather inhomogeneous and far from being complete. 
Most stars, especially in MS, do not show a permanent pe- 
riodic modulation, and therefore deriving a rotation period 
is sometimes a matter of cha nce, unless stars are observed 
for several seasons (see, e.g.. iParihar et al.1 12009). Often a 
beat period is detected and prolonged monitoring is required 
to remove aliasing effects. The effect of large spotted areas 
localised on opposite stellar hemispheres, producing a rota- 
tional modulation with half the true period, should also be 
taken into account (see, e.g.. ICollier Cameron et "al"1 F2009), 

The dispersal of clusters' stars into the field may add 
other selection effects. In order to assemble a complete sam- 
ple one should, in principle, collect data for all stars in each 
cluster as well as data for all cluster's stars that have been 
dispersed into the field, a task which is far from being accom- 
plished. Regarding a Per, which is particularly important to 
constrain the rotational evolution just after the ZAMS, an 
alternative explanation for the presence of a high percentage 
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Figure 7. Cumulative distributions for the equatorial rotational 
velocities derived from the vsini data (thin line) and from the 
rotation period data (thick line) in a Per. 



of fast and ultra-fast rotators is therefore that the dataset 
available may miss information about a significant fraction 
of slower rotators, either because of observational biases or 
because of an intrinsic incompleteness of the sample due to 
stars' dispersal into the field. We note in passing that the 
fraction of ultra-fast rotators (i.e., P ro t < 2 d) in the CDF 
of AB Dor / Pleiades was reduced from ~ 0.2 to the cur- 
rent ~ 0.12 by the data recently acquired bv lHartman et all 
(2010). The same reasoning may also apply to a Per, M50, 
M35, M 34 and M 11. 

Some indication of possible observational biases in the 
alpha Per periods can be deriv ed from a co mparison with 
v sin i measurements. Using the iGaigel (|l993h inversion pro- 
cedures for a r andom orientation of rot ational axes and 
v sin i data from lStauffer et ah! (|l985l . [l989h . we have derived 
the cumulative distribution of equatorial velocities for alpha 
Per stars in our mass r ange. Such procedure s have been ap- 
plied, for instance, by lOueloz et al.l |l998) in an analysis 
of rotational velocities in the Pleiades. We have then com- 
pared such a distribution with the cumulative distribution of 
equatorial velocities derived from the period measurements 
assuming that the radius of all stars in our mass range is 
approximately 1 Rq (Fig. [7J). 

For the slowest rotators, IStauffer et al. I l| 19851 . 1 19891 ) arc 
able to derive only an upper limit of vsini, which for the 
stars in our mass range is ~ 10 km/s. This implies that 
the cumulative distribution of recovered equatorial veloc- 
ities is rather uncertain at the lower end of the velocity 
range. Numerical experiments with several assumptions on 
the vsini distribution below 10 km/s have shown, however, 
that above v aq ~ 15 km/s the derived cumulative distribu- 
tion function is unaffected by any assumption made on the 
vsini distribution below 10 km/s. On the other hand, the 
expected differences in stellar radii in our mass range have 
negligible effects on the cumulative distribution of equato- 
rial velocities derived from the period measurements. There- 
fore, a meaningful comparison can be done between the frac- 
tions of stars with rotational velocities lower than, say, 30 
km/s in both distributions, which are approximately 80% 
according to the v sin i measurements and 50 % according to 
the rotational period measurements. Furthermore, from the 
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IChandrasekhar fc Miinchl l|l95Cj ) relationships (the overline 
denotes average) 



v sm i = —v 
4 



and 



(v sini) 2 



16 



v sini)' 



we derive v — 35.5 km/s, and a 2 = 2800 (km/s) 2 , while the 
rotational velocities inferred from the period measurements 
have v = 63.9 km/s and a 2 = 3700 (km/s) 2 . Finally, the KS 
test applied to these cumulative distribution functions gives 
a probability of only 0.008 that they are drawn from the 
same distribution. We conclude that there is a convincing 
evidence of a bias in the measured period distribution in a 
Per towards an excess of fast rotators. Having carried out 
such a test on the problematic period distribution of a Per, 
a similar comparison for other clusters in our sample would 
be desirable but it is outside the scope of this paper and is 
deferred to future work. 

The DZM with a power-law prescription for r c is found 
quite effective in describing the rotational evolution from 
AB Dor / Pleiades to M 37. Between these two, however, 
it predicts more rotators with P rot > 5 d than observed in 
M 50 and M 35. For M 34 (which seems to have also a higher 
fraction of fast rotators than predicted) and Mil the num- 
ber of available periods seems too low to draw any definitive 
conclusion. Assuming that the DZM contains the essential 
physics to describe the evolution from AB Dor / Pleiades to 
M 37 and that the age estimate is not too poor, one may 
reach the conclusion that a significant fraction of rotators 
with P rot > 5 d are missing in M 50 and M 35. This conclu- 
sion is also supported by a crude comparison between the 
slower rotator tail in the observed CDFs, from which M 50 
and M 35 would appear rotationally older than M 37. At 
ages between the ZAMS and AB Dor /Pleiades, the DZM 
with a power-law prescription for r c is unable to describe the 
high percentage of fast and ultra-fast rotators in a Per, as 
discussed above, and is not able to fit the Tuc-Hor/Col/Car, 
despite the shape of the observed and fitted CDFs' are not 
too different. For these latter associations, however, it should 
be noted that, besides the issue of completeness discussed 
above, the modelling is particular sensitive to age uncertain- 
ties in this evolution phase and to the a Per constraint at 
70 Myr. 

The fit using the two-value prescription for r c is more 
affected by the high percentage of (ultra) fast rotators in 
a Per than the power-law prescription. This prevents from 
obtaining a satisfactory fit for older systems, where such a 
high percentage of fast and ultra-fast rotators is not present. 

Only for M 34, which also seems to show a high per- 
centage of fast rotators, Pks,i is close to unity. Also for Tuc- 
Hor/Col/Car the two-value prescription for r c produces a 
synthetic CDF which is slightly higher than the power-law 
prescription for the faster rotators partly because of the con- 
straint imposed by a Per. 



6 CONCLUSIONS 

We studied the rotational evolution of solar-mass stars from 
PMS to the solar age, using the most complete sample of ob- 



servations on rotation period distributions available to date, 
filling the gap between ~ 6 and 70 Myr with recently ac- 
quired data and implementing a novel prescription on the 
internal rotational coupling in the framework of the DZM. 

Although the completeness and accuracy of observa- 
tional data have considerably improved in recent years, pos- 
sible selection effects (e.g. scarce sensitivity to rotational 
modulation with P rot > 10 d, limits for detection of fast 
or slow rotation periods due to insufficient time sampling 
or extension of the observations), as well as the uncertain- 
ties in clusters' ages, are still severe limitations for the use 
of the period distributions to study stellar rotational evo- 
lution. Nevertheless, we investigated the dependence of the 
core-envelope coupling timescale r c on stellar rotation. 

After the ZAMS, a power-law dependence of r c on the 
stellar rotational velocity produces CDFs whose shapes are 
more similar to observations than the two- value prescription 
and gives an excellent fit for AB Dor / Pleiades and M37. 
While M 34 and Mil have too few period measurements to 
draw firm conclusions, the discrepancies in a Per, M 50, and 
M 35 could be due either to observational biases or to envi- 
ronmental effects, such as those causing early destruction of 
circumstellar discs. 

Comparing the results obtained with the two- value pre- 
scription for t c with those obtained with the power-law, it 
appears that the former better reproduce the extrema while 
the latter the shape of the observed period distributions. 
This conclusion may be a useful guide to formulate more 
sophisticated models of rotational evolution in the future. 

Part of the discrepancy between the models and the 
observations here may also arise from the assumption of a 
single disk lifetime (raise = 5.8 Myr). Assuming instead a 
distribution of disk lifetimes would yield more fast rotators 
at a given age and a broader rotation period distribution at 
the ZAMS, but including this feature in the model is beyond 
the scope of this paper. 

To put our conclusions on a firm ground, however, more 
observations of clusters and associations in the age range 
30 — 100 Myr are required, together with an improvement 
of the model, e.g., by including correlations with the global 
properties of each cluster. 
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